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Abstract. We describe high order accurate and stable fully-discrete finite 
difference schemes for the initial-boundary value problem associated with the 
magnetic induction equations. These equations model the evolution of a mag- 
netic field due to a given velocity field. The finite difference schemes are 
based on Summation by Parts (SBP) operators for spatial derivatives and a 
Simultaneous Approximation Term (SAT) technique for imposing boundary 
conditions. We present various numerical experiments that demonstrate both 
the stability as well as high order of accuracy of the schemes. 



1. Introduction 
In this paper, we study the magnetic induction equation 

(1.1) dtB + curl(B X u) = -udiv(B), 

where the unknown B = B(x, G describes the magnetic field of a plasma in 
three space dimensions with coordinate x = (x^y^z). The above equation models 
the evolution of the magnetic field in the plasma which is moving with a prescribed 
velocity field u(x,t). These equations arise in a wide variety of applications in 
plasma physics, astrophysics and electrical engineering. One important application 
are the equations of magneto- hydro dynamics (MHD). Observe that, by taking 
divergence on both sides of (1.1) we get 

(1.2) (div(B))^ + div (udiv(B)) = 0. 

Hence, if div(Bo(x)) = 0, also div(B(x, t)) = for t > 0. 

There are many forms of induction equations available in literature (see [ , ]). 
Here we are going to work with the following "conservative" symmetric form, 

(1.3) dtB + (u • V) B = M{Du)B, 

where the I^u denotes the gradient of u and the matrix M{Du) is given by 

(—dyU^ — dzU^ dyU^ dzU^ \ 

dxU^ dyU^ —dxU^ — dyU^ j 

We are aware of some other results in the literature related to induction equation 
[1, 2, 3, 4, 11, 6]. But, boundary conditions were not considered either of the 
aforementioned papers. In [ ], authors have described a high order accurate and 
stable finite difference schemes for the initial-boundary problem associated with the 
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magnetic induction equations. The approach is based on a "semi-discrete" approx- 
imation where one discretizes the spatial variable, thereby reducing the equations 
to a system of ordinary differential equations. However, we stress that for nu- 
merical computations also this set of ordinary differential equations will have to 
be discretized in order to be solved. Thus in order to have a completely satisfac- 
tory numerical method, one seeks a fully discrete scheme that reduces the actual 
computation to a solution of a finite set of algebraic equations. 

Our aim in this paper is to design stable and high-order accurate "fully-discrete" 
schemes for initial-boundary value problems corresponding to the magnetic induc- 
tion equations by discretizing the non- conservative symmetric form (1.3). The 
spatial derivatives are approximated by second and fourth-order SBP (Summation- 
By-Parts) operators. The boundary conditions are weakly imposed by using a SAT 
(Simultaneous Approximation Term) and backward Euler method used for tempo- 
ral discretization. The SBP-SAT framework has been used to obtain stable and 
accurate high order schemes for a wide variety of hyperbolic problems in recent 
years. See [ ] and the references therein for more details. 

The rest of this paper is organized as follows: In Section 2, we state the energy 
estimate for the initial-boundary value problem corresponding to (1.3. In Section 3, 
we present the fully-discrete SBP-SAT scheme and show stability. Numerical ex- 
periments are presented in Section 4 and conclusions are drawn in Section 5. 

2. The Continuous problem 

For ease of notation, we shall restrict ourselves to two spatial dimensions in the 
remainder of this paper. Extending the results to three dimensions is straightfor- 
ward. 

In two dimensions, the non- conservative symmetric form (1.3) reads 
(2.1) Bt + AiB^ + A2By - CB = 0, 

where 



with B = {B^,B^) and u = {u^,u^ ) denoting the magnetic and velocity fields 
respectively. In component form, (2.1) becomes 

{B')i + u\B% + u\B% = -{u\B' + (u\B^ 

iB\ + u\B%+u\B% = iu^),B' -{u%B^. 

To begin with, we shall consider (2.1) in the domain (x^y) G 1^ = [0, 1]^. 
We augment (2.1) with initial conditions, 

(2.3) B(x,0)=Bo(x) xen, 

and Dirichlet boundary conditions, 
(2.4) 

l{ni(o,2y,^)>o} (b(0, y, t) = g(0, y, t)^ , l{ni(i,2y,t)<o} (b(1, y, t) = g(l, y, t) 



l{n2(x,0,^)>0} 



{b{x, 0, t) = g(x, 0, t)^ , l{u^x,i,t)<o} (b(x, 1, t) = g(x, 1, t)^ 



where 1a denotes the characteristic function of the set A. Note that we only impose 
boundary conditions on the set where the characteristics are entering the domain. 
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We shall always assume that the initial and boundary data satisfy the compatibility 
conditions, i.e., specific criteria that guarantee smoothness of the solution, see [7]. 

Theorem 2.1. Assume that Bq G H^{n), that g G H^{dn x [0,T]) for T > 
and that and v? are in H^{Q x [0,T]). Then there exists a function B G 
C{[0,T],L'^{n))nL'^{[0,T];H^{n)) which is the unique weak solution of (2.1) with 
the initial and boundary conditions (2.3) and (2.4)- 

Furthermore, it satisfies the following stability estimate 

(2.5) ||B(-,t)||^,(^) 

where a is a positive constant. 

3. FULLY-DISCRETE SCHEME 

To simplify the treatment of the boundary terms we let the computational do- 
main Q be the unit square. It is straightforward to generalize our results to other 
domains by coordinate transformations (see [. ]), and to three dimensions. 

The SBP finite difference schemes for one-dimensional derivative approximations 
are as follows. Let [0, 1] be the domain discretized with Xj = jAx, j = 0, . . . , A^ — 1. 
A scalar grid function is defined as = {wq^ ...wn-i)- To approximate dxW we use a 
summation- by-parts operator Dx = P~^Qx, where Px is a diagonal positive N x N 
matrix, defining an inner product 

1/2 

such that the associated norm ||i^||p^ = {w,w)p^ is equivalent to the norm \\w\\ = 
(Ax k;^)^/^. Furthermore, for Dx to be a summation-by-parts operator we 
require that 

Qx + Qx = Rn — Ln-i 
where Rn and Ln are the N x N matrices: diag (0, . . . , 1) and diag (1, . . . , 0) re- 
spectively. Similarly, we can define a summation- by-parts operator Dy = Py^Qy 
approximating dy. Later we will also need the following Lemma, proven in [^J. 

Lemma 3.1. Let u be a smooth grid function. Then 

(3.1) \\Dx(uow)-uoDxw\\p^ < C'II^:c^IIl-([o,i]) WHp, 

where {u o v)j = UjVj- 

Next, we move on to the two-dimensional case and discretize the unit square 
[0,1]^ using NM uniformly distributed grid points {xi^yj) = {iAx^jAy) for i = 
0, . . . , TV - 1, and j = 0, . . . , M - 1, such that {N - 1) Ax = (M - 1) A^ = 1. We 
order a scalar grid function w{xi^yi) = Wij as a column vector 

W = (^0,0,^0,1, . . . ,^0,(M-1),^1,0, ...... . ,^(Ar-l),(M-l))^ • 

To obtain a compact notation for partial derivatives of a grid function, we use 
Kronecker products. The Kronecker product of an A^i x 7V2 matrix A and an 
Ml X M2 matrix B is defined as the A^iMi x N2M2 matrix 

/ aiiB . . . aiN^B 

(3.2) A®B= : : 

\aN^iB ... aN-^N^Bj 
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For appropriate matrices A, C and the Kronecker product obeys the fohowing 
rules: 

(3.3) {A (8) B){C ^D) = {AC (g) BD), 

(3.4) (A (8) 5) + (C (g) L>) = (A + C) (B + D), 

(3.5) {A^Bf = {A^ 0B^). 

Using Kronecker products, we can define 2-D difference operators. Let In denote 
the n X n identity matrix, and define 

= Im, = In ^ Dy. 

For a smooth function w{x,y), {'Oxw)ij ^ dxw{xi,yj) and similarly {dyw)ij ^ 
dyw{xi,yj). 

Set V = Px® Py, define {w,v)r = w^Fv and the corresponding norm \\w\\p = 

1/2 

{w, w)-p . Also define IZ = Rn^Im^ ^ = Ln^Im^ U = In^Rm and P = I^^Lm- 
For a vector valued grid function V = (V-^, y^), we use the following notation 



dxy = 



and so on. In the same spirit, the P inner product of vector valued grid functions is 
defined by (V, W)p = {V^,W^)r + (V^"^2)p. We also introduce (a small) time 
step At > 0, and use the notation 

for any function p : [0, T] R. Write = nAt forn G No = N U {0}. We wiU use 
the notation V^{xi,yj,t^) = V^j'^ and so on. 

Remark 3.1. Note that the Kronecker products is just a tool to facilitate the no- 
tation. In the implementation of schemes using the operators in the Kronecker 
products we can think of these as operating in their own dimension, i.e., on a spe- 
cific index. Thus, to compute dxW, we can view w as a field with two indices, and 
the one- dimensional operator Dx will operate on the first index since it appears in 
the first position in the Kronecker product. 

The usefulness of summation by parts operators comes from this lemma. 

Lemma 3.2. For any grid functions v and w, we have 



{v, dxw)p + {dxV, w)p = [(7^ - C) {In ^ Py)] w 

•6) 

{v,dyw)p + {dyV,w)p = V [{U -V){Px (8)/m)]'^- 
Observe that this lemma is the discrete version of the equality 

// v{dxw)dxdy^ // {dxv)wdxdy= / v{l,y)w{l,y) - v{0,y)w{0,y) dy. 
J Jn J Jn Jo 



Proof. We calculate 



{V, dxW) = V' {Px Py) [Px Qx ^ Im) 



w 



= v^Qx (8) PyW 

= -v^Ql (8) PyW + v^{Qx + Ql) (8) PyW 

= -{Px^Qx (8 iMvf {Px ^Py)w^ {Rn - Ln) (8 PyW 
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= - i^xvf {Pa: ^Py)w^ (7^ - C) {In ^ Py)w. 

The second equality is proved similarly. □ 

Before we define our numerical schemes, we collect some useful results in a 
lemma. 

Lemma 3.3. If u is a grid function, then 

(V, u o c),V)p = \\^[{n - C){In Py)] {u o V) 

+ i(iioc)^V-c)^ (ixoV),V)p, 

(3.7) 2 

(V, u o c)^V)p = i^^^m - ^){Px ^ Im)] {u o V) 

Proof. To show (3.7), first note that since P is diagonal, {u o c)V, V)p = (c)V,i^ o 
V)p. We use Lemma 3.2 to calculate 

{u o c)^V, V)p = (c)^(^/ o V) , V)p + {u o c)^V - c)^ (^1 o V) , V)p 

= -{uo\^ c), V)p + (7^ - £) {In ^ Py) {u o V) 

+ (ix o c)^ V - (iioV),V)p 
= - (ii o c)^ V, V)p + (7^ - £) (/tv ^ Py) {u o V) 

+ (^x o c)^, V - c)^ (x/oV),V)p. 

This shows the first equation in (3.7), the second is proved similarly. □ 

Now we are in a position to state our scheme(s). For £ = 1 or 2 we will use 
the notation for both the grid function defined by the function u^{x^y) and for 
the function itself. Similarly, for the boundary values, we use the notation h and g 
for both discrete and continuously defined functions. Hopefully, it will be apparent 
from the context what we refer to. 

The differential equation (2.1) will be discretized in an obvious manner. We 
incorporate the boundary conditions by penalizing boundary values away from the 
desired ones with a (9(1/ Ax) term. To this end set 

B = [{P-^ Im) {^cC + E7^7^) + {In P"') {^vV + ^uU)] , 

where S/:, T^^^ T^t> and T^u are diagonal matrices, with components {(Jc)jj ordered 
in the same way as in ((3.2)) (and similarly for the other penalty matrices), to be 
specified later. 

With this notation the scheme for the differential equation (2.1) reads 

(3.8) L)+V^ + xx^'^+^oc)^V^+^+ii2'^+^oc)^V^+^-C^+^V^+^ =S(V^+^-g^+^), 
while is given. Here denotes the matrix 
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Theorem 3.1. Let V be as solution to (3.8) with g = 0. // the constants in B is 

chosen as 

(3.9) 

[crn)N-ij < , [crc)oj < ^ , [cru)i,M-i < , 

.... ^^'+(^i,0) 
and [aT))i,o < ^ , 

then 

(3.10) l|V"||^ <e^^||V°||p, 

where v}'^ = {u^ V 0)^ u^'~ = {u^ A 0)^ for / = 1, 2. K is a constant chosen in such 
a way that (1 — cAt)~-^ < (1 + i^At) for sufficiently small At, where c is a constant 
depending on , , and their derivative approximations, hut not on N or M. 

Proof Taking the P inner product of (3.8) and V"^+^, we get 

I l|V"+ip - I l|V"||^ + ^ ||V"+i - V"||p = -At (V"+\ui'»+i oO,V"+i)p 
- At (V"+\ o OyV"+i)p + At (v"+\ C"+iV"+i)p + At (V"+\ SV"+i)p , 

Using Lemma 3.3 we get 



l|V"||^ + -||V»+i-V"| 
= -At^{V"+Y [{TZ - C){In ® Py)] o V"+i) 

- At^iV^+Y W - ^){Px ® Im)] o V"+i) 

- At^ oO^V"+^ o V"+i),V"+i)p 

- At^ oj)j,V"+i -0j,(u2'"+^ o V"+i),V"+i)p 

+ At (V"+\ CV"+^)p + At (V"+\BV"+i)p . 

Note that by (3.1), 

I o c),V"+^ - c)^(zi^'"+i o V"+^), V"+^)p| < c II V"+^ lip , 
(3.11) |(w2.«+io0j,V"+^ -0j,(u2'"+ioV"+^),V"+^)p| <c||V"+^||p, 

|(V"+\CV"+i)p| < c||V"+i||p, 

for some constant c depending on the first derivatives of v} and v? . Using the 
conditions (3.9) we arrive at 

||V"+i||p < ||V"||^+cAt||V"+i||p 

Now we can use the fact that (1 — cAt)"-*- < (1 + if At) for sufficiently smah At. 
Consequently this gives the required bound (3.10). □ 
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4. Numerical Experiment 

We test the fully-discrete SBP-SAT scheme of the previous section on a suite of 
numerical experiments in order to demonstrate the effectiveness of these schemes. 
We will use two different schemes : SBP2 and SB PA scheme which are second-order 
(first-order) and fourth order (second-order) accurate in the interior (boundary) 
resulting in an overall second and third-order of accuracy. 

In this experiment, we consider (2.1) with the divergence- free velocity field 
u{x^y) = {—y^x)^. The exact solution can be easily calculated by the method 
of characteristics and takes the form 

(4.1) B(x,t) = i?(t)Bo(i?(-t)x), 

where R{t) is a rotation matrix with angle t and represents rotation of the initial 
data about the origin. 

We consider the same test setup as in [ ] and [ ] by choosing the divergence 
free initial data, 

(4.2) Bo(a:,y)=4(^-^)e-2°((-V2)W), 

and the computational domain [—1,1] x [—1, 1]. Since the exact solution is known 
in this case, one can in principle use this to specify the boundary data g. Instead, 
we decided to mimic a free space boundary (artificial boundary) by taking ^ = 0. 
(which is a good guess at a far- field boundary). 

We run this test case with SBP2 and SBPA schemes and present different sets 
of results. In Figure 4.1, we plot |B| = (\B^\^ + |52|2^)i/2 ^-^^g ^ ^ ^ ^Yv^li- 
rotation) and t = 27r (one full rotation) with the SBP2 and SBPA schemes. As 
shown in this figure, SBP2 and SBPA schemes resolve the solution quite well. In 
fact, SBPA is very accurate and keeps the hump intact throughout the rotation. In 



Grid size 


SBP2 


rate 


SBPA 


rate 


40x40 


6.9 xlO^ 




8.0 xlO^ 




80x80 


2.1 xlO^ 


1.7 


5.0 xlO-^ 


4.0 


160x160 


5.5 xlO^ 


2.0 


4.5 xlO-2 


3.5 


320x320 


1.3 xlO° 


2.0 


5.1 xlO-^ 


3.1 


640x640 


3.3 xlO-i 


2.0 


6.4 xlO-"^ 


3.0 



Table 4.1. Relative percentage errors in for |B| at time t = 2tt 
and rates of convergence with SBP2 and SB PA schemes. 

Table 4.1, we present percentage relative errors in P. The errors are computed at 
time t = 27T (one rotation) on a sequence of meshes for both the SBP2 and SBPA 
schemes. The results show that the errors are quite low, particularly for SBPA 
and the rate of convergence approaches the expected values of 2 for SBP2 and 3 
for SBPA. Furthermore, the order of accuracy is unaffected at these resolutions by 
using zero Dirichlet boundary data instead of the exact solution at the boundary. 

5. Conclusion 

We have considered a fully-discrete scheme for the magnetic induction equations 
that arise as a submodel in the MHD equations of plasma physics. In future, our 
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I 



(a) half rotation, SBP2 



I 



(b) full rotation, SBP2 



(c) half rotation, SBP4 



(d) full rotation, SBP4 



Figure 4.1. Numerical results for IBI. 



plan is to extend the semi-discrete scheme given in [o] to a semi-implicit fully- 
discrete scheme. We would like to show the stabilty of the aforementioned semi- 
implicit scheme in case of magnetic induction equations with resistivity. 
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